#Function to plot the marginal effects and the differences
marginal_effect_plot <- function(dataframe, facet = TRUE, models = FALSE, ORHR = FALSE, base_size = 16) {
if (is.null(dataframe$name)) {
dataframe <- dataframe %>%
pivot_longer(!!names(.[3]), values_drop_na = T)
}
if (facet == TRUE & models == FALSE & ORHR == FALSE) {
plot <- ggplot(dataframe, aes(estimate, str_replace_all(term, output_names), xmin = conf.low, xmax = conf.high,
col = factor(value), linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
group = value))+
geom_linerange(position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 2) +
geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
facet_wrap(~str_replace_all(name, output_names))+
labs(x = "Marginal Effect", y = NULL)+
theme_base(base_size = base_size)+
scale_color_viridis(discrete = T, option = "D", name = NULL)+
scale_linewidth_discrete(range = c(0.4, 1.2))+
scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
guides(color = guide_legend(reverse=TRUE),
linewidth = "none")+
theme(legend.position = "top",
legend.box = "vertical", plot.background=element_blank())
print(plot)
} else if (facet == TRUE & models == TRUE & ORHR == FALSE) {
plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high,
col = model, linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
group = model))+
geom_linerange(position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 2) +
geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
facet_wrap(~str_replace_all(name, output_names))+
labs(x = "Marginal Effect", y = NULL)+
theme_base(base_size = base_size)+
scale_color_viridis(discrete = T, option = "D", name = NULL)+
scale_linewidth_discrete(range = c(0.4, 1.2))+
scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
guides(color = guide_legend(reverse=TRUE),
linewidth = "none")+
theme(legend.position = "top",
legend.box = "vertical", plot.background=element_blank())
print(plot)
} else if (facet == TRUE & models == TRUE & ORHR == TRUE) {
plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high,
col = model, linewidth = factor(c_lvl, levels = c("0.95", "0.9")),
group = model))+
geom_linerange(position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 2) +
geom_vline(xintercept = 1, color = "darkred", linetype = 2)+
facet_wrap(~str_replace_all(name, output_names))+
labs(x = "Marginal Effect (Odds/Hazard Ratio)", y = NULL)+
theme_base(base_size = base_size)+
scale_color_viridis(discrete = T, option = "D", name = NULL)+
scale_linewidth_discrete(range = c(1.2, 2.2))+
scale_x_continuous(breaks = seq(0,6, by=0.5))+
guides(color = guide_legend(reverse=TRUE),
linewidth = "none")+
theme(legend.position = "top",
legend.box = "vertical", plot.background=element_blank())
print(plot)
} else {
plot <- ggplot(dataframe, aes(estimate, factor(value), xmin = conf.low, xmax = conf.high,
group = factor(value),
linewidth = factor(c_lvl, levels = c("0.95", "0.9"))))+
geom_linerange(position = position_dodge(0.4)) +
geom_point(position = position_dodge(0.4), size = 3) +
geom_vline(xintercept = 0, color = "darkred", linetype = 2)+
labs(x = "Marginal Effect", y = NULL)+
theme_base(base_size = base_size)+
scale_color_viridis(discrete = T, option = "D", name = NULL)+
scale_linewidth_discrete(range = c(1.4, 2.4))+
scale_x_continuous(breaks = seq(-0.2,0.2, by=0.05))+
guides(color = guide_legend(reverse=TRUE),
linewidth = "none")+
theme(legend.position = "top",
legend.box = "vertical", plot.background=element_blank())
print(plot)
}
}
#### Lists of names ####
output_names <- c("GOV_turnover_lead" = "Change of government (t+1)",
"GOV_turnover_lagged" = "Change of government (t-1)",
"GOV_turnover" = "Change of government",
"minister_turnover" = "Change of minister",
"bureaucratic_buffer_i" = "Bureaucratic buffer & position type",
"bureaucratic_buffer" = "Bureaucratic buffer",
"organizational_buffer2" = "Organizational buffer",
"state_secretaries_NSDID" = "Political buffer",
"state_secretaries" = "Political buffer",
"years_since_last_turnover_2" = "Time since change of government",
"years_since_last_turnover_2_categories" = "Time since change of government categorial",
"match_start_end_PMparty_flipped" = "Party incongruence",
"election_year" = "Election year",
"age_year" = "Age",
"gender" = "Gender",
"temporary_employment" = "Temporary appointment",
"pre_unionsoppløsning" = "Before-1906",
"foreign_affairs" = "Ministry of Foreign Affairs",
"government_PMparty_share_posts_end_of_year" = "PM-party share of cabinet",
"years_since_DEP_start" = "Time since ministry created",
"dep_end" = "Ministry terminated",
"pol_flexible" = "Exposed to change of government",
"NSD_ID" = "Ministry Fixed Effects",
"decade" = "Decade",
"decade2" = "Decade",
"duration" = "Senior bureaucrat tenure in years",
"WWII" = "WWII",
"age_group_year2" = "Age Group",
"multiple_L1" = "Multiple Level-1 bureuacrats",
"dep.rad_turnover" = "Senior bureaucrat turnover")
moderators_names <- str_replace_all(sapply(c("Baseline",
"government_PMparty_share_posts_end_of_year",
"years_since_DEP_start", "dep_end",
"multiple_L1", "age_group_year2"),
function(x){rep(x, 6)}), output_names)
DV_names <- lapply(c("Senior bureaucrat turnover",
"Including ministry switching",
"Including ministry switching \n and promotion/demotion"),
function(x){rep(x, 6)}) %>% unlist()
#### modelsummary fit statistics ####
gm <- tribble(~raw, ~clean, ~fmt,
"FE: time", "FE: Time", 0,
"FE: NSD_ID", "FE: Ministry", 0,
"FE: decade", "FE: Decade", 0,
"nobs", "N", 0,
"r.squared", "R2", 2,
"aic", "AIC", 2)
###################################
####preparing data for analysis####
###################################
#Load data
df <- readRDS(file = "data/TCS_NOR.RDS")
#recoding
df$match_start_end_PMparty_flipped <- if_else(df$match_start_end_PMparty == 1, 0, 1)
df$election_year <- if_else(is.na(df$election_year), 0, df$election_year)
df$military <- ifelse(is.na(df$military), 0, df$military)
#### creating subsets of the data ####
raw <- df %>% filter(year >= 1884, military == 0) #removing military DGs and observations before parliamentarism
#pensioners are removed | and WWII
df <- raw %>%
mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>%
filter(age_year <= 70, WWII == 0)
#pensioners are removed
df_med_WWII <- raw %>%
mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>%
filter(age_year <= 70)
#spells must be more than 2 years and pensioners are removed
df2 <- df %>%
filter(less_than_two_years == 0) %>%
mutate(dep.rad_turnover = if_else(pension_age_year == 1, 0, dep.rad_turnover)) %>%
filter(age_year <= 70)
#only Level-1 bureaucrats
df_level1 <- df %>% filter(bureaucratic_buffer == 0)
#### dataframe with all used variables
sub_df <- df_med_WWII %>%
select(start, stop, dep.rad_turnover, GOV_turnover, years_since_last_turnover_2, match_start_end_PMparty_flipped,
election_year, age_year, gender, temporary_employment, pol_flexible, pre_unionsoppløsning, WWII, decade2,
dep_end, years_since_DEP_start, government_PMparty_share_posts_end_of_year,
foreign_affairs, WWI, post_WWII, NSD_ID, bureaucratic_buffer, year,
state_secretaries, proximity_to_GOV_turnover,proximity_to_GOV_turnover_2_categories, years_since_last_turnover_2_categories,
GOV_turnover_any_change, GOV_turnover_new_PM, GOV_turnover_new_Parties,
GOV_turnover_during_term, GOV_turnover_election, government_percent_parliament_end_of_year,
government_single_party_end_of_year, government_minority_end_of_year, government_left_right_gov_end_of_year,
gov_turnover_during_spell, pension_age_propensity_age_year, pension_age_propensity_age_at_start,
match_start_end_party, match_start_end_left_right, match_start_end_coalition_party,
minister_turnover, political_career,
GOV_turnover_lagged, GOV_turnover_lead, id_spells, duration, bureaucratic_buffer_i,
y_since_PS_in_min, decade_since_PS_in_min,rounded_y_since_PS_in_min, organizational_buffer,
organizational_buffer2, after_political_appointees_NSDID, after_state_secretaries_NSDID, after_poladv_NSDID,
next_position_sector, education_main, education_level, post_SLS, pension_age_propensity_age_year,
age_group_year2, multiple_L1, ministry_portfolio)
####Figure 9 - Types of government in Norway 1884-2024 ####
governments <- readRDS("data/Governments_NOR.RDS") %>%
mutate(government_type = case_when(regjeringstype == "Flertall ettparti" ~ "Single-Party Majority",
regjeringstype == "Flertall minste vinnende koalisjon" |
regjeringstype == "Flertall overtallig koalisjon"  ~ "Multi-Party Majority",
regjeringstype == "Mindretall ettparti" ~ "Single-Party Minority",
regjeringstype == "Mindretall koalisjon" ~ "Multi-Party Minority"))
#time labels
table <- governments %>%
group_by(government_type) %>%
summarise(time_with_government_type = sum(lengde_nesten_alle_endringer, na.rm = F),
time_with_government_type_percentage = time_with_government_type /
(as.numeric( ymd("2024-12-31") - ymd("1884-06-26") ) ))
time_labels <- paste(table$government_type,
paste0( paste0("(", scales::percent( as.numeric( table$time_with_government_type_percentage) ) ), ")" )
, sep = " ")
#plot
obj <- ggplot(governments,
aes(xmin = tiltrådt,
xmax = avgang,
ymin = 0,
ymax = percent,
fill = government_type))+
geom_rect()+
geom_hline(yintercept =0.5, linetype = 2, col = "red")+
scale_x_date(expand = c(0.025, 0), date_breaks = "5 years", date_labels = "%Y", limits = ymd(c("1883-01-01", "2025-01-31") ) )+
scale_y_continuous(limits = c(0, 1), breaks = seq(0, 1, by = 0.10), labels = scales::percent)+
scale_fill_viridis(discrete = T, labels = time_labels)+
theme_base(base_size = 16)+
labs(fill = NULL, y = "Percentage of seats in parliament", x = NULL) +
theme(legend.position = "bottom",
legend.box = "vertical",
plot.background=element_blank(),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))+
guides(fill=guide_legend(nrow = 3,
byrow = T,
title.position = "top"
))
print_plot(plot = obj, name = "governmentsovertime")
##Time with different governments
#minority
18.5+36.4
#single-party
23.8+36.4
##number of change in prime minister and party
sum(governments$statsminister_parti != lag(governments$statsminister_parti), na.rm = T)
# number of change after an election
sum(governments$lag_årsak_til_avgang[governments$statsminister_parti != lag(governments$statsminister_parti)] != "Valg", na.rm = T)
# Type of governments formed
table(governments$minority_majority[governments$statsminister_parti != lag(governments$statsminister_parti)])
# Type of governments formed
table(governments$koalisjon_ettparti[governments$statsminister_parti != lag(governments$statsminister_parti)])
####Figure 10 - The Percentage of ministries employing at least 1 state secretary and/or 1 permanent secretary over time ####
plotdata <- df %>%
distinct(year, NSD_ID, PS_in_dep, state_secretaries_NSDID, GOV_turnover) %>%
reframe(share_dep_rad = sum(PS_in_dep == 1)/n(),
share_statesec = sum(state_secretaries_NSDID == 1)/n(),
GOV_turnover_year = if_else(GOV_turnover == 1, year, NA),
.by = year) %>%
pivot_longer(cols = share_dep_rad:share_statesec)
variance_in_positions <- ggplot(plotdata, aes(year, value, col = name))+
geom_segment(aes(x = GOV_turnover_year, y = 0, yend = 0.03), col="black")+
geom_line(size=1.3)+
ylab("Percentage of Ministries with Position") +
xlab(NULL) +
theme_base(base_size = 26)+
scale_color_viridis(discrete = T, option = "D", name = NULL, labels = c("Permanent Secretary", "State Secretary"))+
scale_y_continuous(labels = scales::percent)+
scale_x_continuous(breaks = seq(1884,2025, by=20))+
guides(color = guide_legend(reverse=TRUE))+
theme(legend.position = "top",
legend.box = "vertical", plot.background=element_blank())
print_plot(variance_in_positions, name = "variance_in_positions")
####Table 3-5 - Descriptive Statistics ####
descriptive_stats <- sub_df %>%
group_by(id_spells) %>%
dplyr::select(duration, dep.rad_turnover, GOV_turnover, bureaucratic_buffer, state_secretaries,
election_year, age_year, gender, temporary_employment, pol_flexible, foreign_affairs,
pre_unionsoppløsning, government_PMparty_share_posts_end_of_year,
years_since_DEP_start, dep_end) %>%
mutate(duration = round(as.numeric(duration/365.25))) %>%
ungroup() %>%
dplyr::select(-id_spells)
names(descriptive_stats) <- str_replace_all(names(descriptive_stats), output_names)
#All
stargazer(as.data.frame(descriptive_stats),
omit.summary.stat = c("p25", "p75"),
title = "Descriptiv statistics: All senior bureaucrats",
align=F,
no.space=T,
model.numbers=T,
label = c("Tab:desc1"),
type = "latex",
out = output(name = "descriptivstats_complete", format = "tex")
)
##Level-1
highdata <- as.data.frame(descriptive_stats %>%
filter(`Bureaucratic buffer` == 1) %>%
select(-`Bureaucratic buffer`))
stargazer(highdata,
omit.summary.stat = c("p25", "p75"),
title = "Descriptiv statistics: Level-1 Bureaucrats",
align=F,
no.space=T,
model.numbers=T,
label = c("Tab:desc2"),
type = "latex",
out = output(name = "descriptivstats_highest", format = "tex")
)
##Level-2
lowdata <- as.data.frame(descriptive_stats %>%
filter(`Bureaucratic buffer` == 0) %>%
select(-`Bureaucratic buffer`) )
stargazer(lowdata,
omit.summary.stat = c("p25", "p75"),
title = "Descriptiv statistics: Level-2 Bureaucrats",
align=F,
no.space=T,
model.numbers=T,
label = c("Tab:desc3"),
type = "latex",
out = output(name = "descriptivstats_lowest", format = "tex")
)
####Figure 11 - The number of senior bureaucrats that exited their position at different intervals of tenure in years####
df$censoring <- if_else(df$Sensurert == "Høyre", "Right-censoring", "Exit")
df_med_WWII$censoring <- if_else(df_med_WWII$Sensurert == "Høyre", "Right-censoring", "Exit")
length <- ggplot(df_med_WWII %>%
distinct(id_spells, duration, censoring),
aes(duration/365.25, fill = censoring))+
geom_bar(width = 0.9)+
scale_x_continuous(limits = c(0, 41), breaks = c(0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40))+
scale_y_continuous(limits = c(0,120))+
theme_base(base_size = 16)+
scale_fill_viridis(discrete = T, name = "")+
labs(x = "Senior Bureaucrat Tenure in Years", y = "Number of Senior Bureaucrats") +
theme(legend.position = "bottom",
legend.box = "vertical", plot.background=element_blank())
print_plot(length, name = "length")
##length of spells at different points in time
df %>%
distinct(id_spells, duration, Sensurert, decade_start) %>%
filter(Sensurert == "ikke sensurert") %>%
mutate(periods = case_when(decade_start < 1940 ~ "Before 1940",
decade_start >= 1950 & decade_start < 1980 ~ "1950-1980",
decade_start >= 1980 & decade_start < 2010 ~ "1980-2010",
TRUE ~ "After 2010")) %>%
summarise(median = median(duration)/365.25,
.by = c(periods))
####Figure 12 - Kaplan-Meier survival curves for senior bureaucrats in years with vs. without change of government, subset by bureaucratic and political buffering####
#Bureaucratic Buffering
km1 <- survfit(Surv(start, stop, dep.rad_turnover) ~
(GOV_turnover) + (bureaucratic_buffer), df, robust = T, id = id_spells)
plotdata1 <- ggsurvplot(km1)$data.survplot
plotdata1$start <- plotdata1$time-1
plotdata1 <- plotdata1 %>%
tidyr::pivot_longer(cols = c(time, start), values_to = "time", names_to = NULL) %>%
mutate(x_variabel = if_else(bureaucratic_buffer == 1, "Level-2", "Level-1"),
grp_variabel = "Bureaucratic buffer") %>%
select(-bureaucratic_buffer)
#Political Buffering
km2 <- survfit(Surv(start, stop, dep.rad_turnover) ~
(GOV_turnover) + (state_secretaries_NSDID), df %>% filter(bureaucratic_buffer == 0), robust = T, id = id_spells)
plotdata2 <- ggsurvplot(km2)$data.survplot
plotdata2$start <- plotdata2$time-1
plotdata2 <- plotdata2 %>%
tidyr::pivot_longer(cols = c(time, start), values_to = "time", names_to = NULL) %>%
mutate(x_variabel = if_else(state_secretaries_NSDID == 1, "Level-1 & Political buffer", "Level-1 & No buffer"),
grp_variabel = "Political buffer") %>%
select(-state_secretaries_NSDID)
#Combining data
plotdata <- full_join(plotdata1, plotdata2)
#Ploting
obj <- ggplot(plotdata, aes(time, surv,
ymax = upper, ymin = lower,
fill = GOV_turnover, col = GOV_turnover))+
geom_line(linetype = 1, linewidth = 1)+
ggh4x::facet_nested_wrap(vars(grp_variabel, x_variabel), nrow = 2, ncol = 2)+
scale_x_continuous(limits = c(0, 41), breaks = seq(0, 40, by = 5))+
scale_y_continuous(limits = c(0,1), labels = scales::percent)+
theme_base(base_size = 16)+
scale_color_viridis(discrete = T, name = "", labels = c("Within-term", "Change of government"), direction = -1)+
scale_fill_viridis(discrete = T, name = "", labels = c("Within-term", "Change of government"), direction = -1)+
labs(x = "Senior Bureaucrat Tenure in Years", y = "") +
theme(legend.position = "bottom",
legend.box = "vertical", plot.background=element_blank())
print_plot(obj, name = "kapmeierplots")
#### Figure 13 - Hazard/Odds ratios with 95% (small stroke width) and 90% (big stroke width) confidence bands for the marginal effect of change of government when estimating the effect with alternative estimators.  ####
## Bureaucratic buffer
#Cox
cox_no_FE <- coxph(Surv(start, stop, dep.rad_turnover) ~ GOV_turnover * bureaucratic_buffer +
election_year + age_year + gender + temporary_employment + pre_unionsoppløsning + as.factor(decade) + as.factor(NSD_ID),
cluster = PersonID,
df)
modelsummary(cox_no_FE, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("cox_no_FE", "html"),
coef_rename = output_names,
coef_omit = "decade|NSD_ID",
exponentiate = TRUE,
gof_map = gm,
add_rows = data.frame(rbind(cbind("FE: Time", " "),
cbind("FE: Ministry", "X"),
cbind("FE: Decade", "X"),
cbind("Number of senior bureaucrats", length(unique(df$id_spells[is.finite(df$bureaucratic_buffer)]))),
cbind("Number of events", cox_no_FE$nevent))),
title = "Cox regression results: DV = Turnover of senior bureaucrat",
notes = list("Note: Hazard Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
#Logistic
GLM_no_FE <- feglm(dep.rad_turnover ~ GOV_turnover * bureaucratic_buffer + time + time_squared + time_cubed +
csw0(election_year, age_year, gender, temporary_employment, pre_unionsoppløsning) | decade + NSD_ID,
cluster = "PersonID",
family = "binomial",
df)
modelsummary(GLM_no_FE, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("GLM_no_FE", "html"),
coef_rename = output_names,
exponentiate = TRUE,
col.names = NULL,
gof_map = gm,
title = "Step-wise Logistic regression results: DV = Turnover of senior bureaucrat",
notes = list("Note: Odds Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
## Political buffer
cox_no_FE_polbuff <- coxph(Surv(start, stop, dep.rad_turnover) ~ GOV_turnover * state_secretaries_NSDID +
election_year + age_year + gender + temporary_employment + pre_unionsoppløsning + as.factor(decade) + as.factor(NSD_ID),
cluster = PersonID,
df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
modelsummary(cox_no_FE_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("cox_no_FE_polbuff", "html"),
coef_rename = output_names,
coef_omit = "decade|NSD_ID",
exponentiate = TRUE,
gof_map = gm,
add_rows = data.frame(rbind(cbind("FE: Time", " "),
cbind("FE: Ministry", "X"),
cbind("FE: Decade", "X"),
cbind("Number of senior bureaucrats", length(unique(df$id_spells[is.finite(df$bureaucratic_buffer) & df$bureaucratic_buffer_i != "Level-2 (Director General)"]))),
cbind("Number of events", cox_no_FE_polbuff$nevent))),
title = "Cox regression results: DV = Turnover of senior bureaucrat",
notes = list("Note: Hazard Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
GLM_no_FE_polbuff <- feglm(dep.rad_turnover ~ GOV_turnover * state_secretaries_NSDID + time + time_squared + time_cubed +
csw0(election_year, age_year, gender, temporary_employment, pre_unionsoppløsning) | decade + NSD_ID,
cluster = "PersonID",
family = "binomial",
df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
modelsummary(GLM_no_FE_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("GLM_no_FE_polbuff", "html"),
coef_rename = output_names,
exponentiate = TRUE,
col.names = NULL,
gof_map = gm,
title = "Step-wise Logistic regression results: DV = Turnover of senior bureaucrat",
notes = list("Note: Odds Ratios with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
plotdata <- bind_rows(
cox_glm_marginal_effects_function(cox_no_FE, "bureaucratic_buffer", mod = "Cox"),
cox_glm_marginal_effects_function(cox_no_FE_polbuff, "state_secretaries_NSDID", mod = "Cox"),
cox_glm_marginal_effects_function(GLM_no_FE[[6]], "bureaucratic_buffer", mod = "Logistic"),
cox_glm_marginal_effects_function(GLM_no_FE_polbuff[[6]], "state_secretaries_NSDID", mod = "Logistic"))
print_plot(marginal_effect_plot(plotdata, models = T, ORHR = T), name = "Coefplot_logistic_cox")
#### Figure 14 - Marginal effects with 95\% (small stroke width) and 90\% (big stroke width) confidence bands for the marginal effect of change of government on alternative specifications of the dependent variable: senior bureaucrat turnover ####
## Bureaucratic buffer
OLS_alternative_DVs <- feols(c(dep.rad_turnover, dep.rad_turnover2, dep.rad_turnover3) ~ GOV_turnover * bureaucratic_buffer +
election_year + age_year + gender + temporary_employment + pre_unionsoppløsning | time+NSD_ID+decade,
cluster = "PersonID", df)
names(OLS_alternative_DVs) <- c("Senior bureaucrat turnover", "Including ministry switching in senior bureaucrat turnover",
"Including ministry switching \n and promotion/demotion in senior bureaucrat turnover")
modelsummary(OLS_alternative_DVs, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("OLS_alternative_DVs", "html"),
coef_rename = output_names,
gof_map = gm,
title = "OLS regression results: Different DV operationalizations",
notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
## Political buffer
OLS_alternative_DVs_polbuff <- feols(c(dep.rad_turnover, dep.rad_turnover2, dep.rad_turnover3) ~ GOV_turnover * state_secretaries_NSDID +
election_year + age_year + gender + temporary_employment + pre_unionsoppløsning | time+NSD_ID+decade,
cluster = "PersonID", df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
names(OLS_alternative_DVs_polbuff) <- c("Senior bureaucrat turnover", "Including ministry switching in senior bureaucrat turnover",
"Including ministry switching \n and promotion/demotion in senior bureaucrat turnover")
modelsummary(OLS_alternative_DVs_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("OLS_alternative_DVs_polbuff", "html"),
coef_rename = output_names,
gof_map = gm,
title = "OLS regression results: Different DV operationalizations",
notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
## Calculating marginal effects and diff-in-diff
Coefplot_OLS_alternative_DVs <- model_marginal_effects(OLS_alternative_DVs) %>% mutate(model = DV_names)
Coefplot_OLS_alternative_DVs_polbuff <- model_marginal_effects(OLS_alternative_DVs_polbuff) %>% mutate(model = DV_names)
## Outputting Coefficient Plots
plotdata <- bind_rows(
Coefplot_OLS_alternative_DVs,
Coefplot_OLS_alternative_DVs_polbuff
) %>%
pivot_longer(c(state_secretaries_NSDID, bureaucratic_buffer), values_drop_na = T) %>%
mutate(model = factor(model, levels = rev(c("Senior bureaucrat turnover",
"Including ministry switching in senior bureaucrat turnover",
"Including ministry switching \n and promotion/demotion in senior bureaucrat turnover"))))
print_plot(marginal_effect_plot(plotdata, models = T, facet = T), name = "coefplot_altDVs")
#### Figure 15 - Marginal effects Additional Moderator Variables ####
## Bureaucratic buffer
OLS_robust_moderators <- feols(dep.rad_turnover ~ GOV_turnover * bureaucratic_buffer +
election_year + age_year + gender + temporary_employment + pre_unionsoppløsning +
sw0(government_PMparty_share_posts_end_of_year, years_since_DEP_start, dep_end, multiple_L1, age_group_year2) | time+NSD_ID+decade,
cluster = "PersonID", df)
names(OLS_robust_moderators) <- str_replace_all(c("Baseline", "government_PMparty_share_posts_end_of_year", "years_since_DEP_start", "dep_end", "multiple_L1", "age_group_year2"), output_names)
modelsummary(OLS_robust_moderators, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("OLS_robust_moderators", "html"),
coef_rename = output_names,
gof_map = gm,
title = "OLS regression results: DV = Turnover of senior bureaucrat",
notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
## Political buffer
OLS_robust_moderators_polbuff <- feols(dep.rad_turnover ~ GOV_turnover * state_secretaries_NSDID +
election_year + age_year + gender + temporary_employment + pre_unionsoppløsning +
sw0(government_PMparty_share_posts_end_of_year,years_since_DEP_start, dep_end, multiple_L1, age_group_year2) | time+NSD_ID+decade,
cluster = "PersonID", df %>% filter(bureaucratic_buffer_i != "Level-2 (Director General)"))
names(OLS_robust_moderators_polbuff) <- str_replace_all(c("Baseline", "government_PMparty_share_posts_end_of_year", "years_since_DEP_start", "dep_end", "multiple_L1", "age_group_year2"), output_names)
modelsummary(OLS_robust_moderators_polbuff, conf_level = 0.95, estimate = "{estimate}", statistic = "[{conf.low}, {conf.high}] {stars}",
output = output("OLS_robust_moderators_polbuff", "html"),
coef_rename = output_names,
gof_map = gm,
title = "OLS regression results: DV = Turnover of senior bureaucrat",
notes = list("Note: Estimate with 95\\% CI (calculated with robust standard errors clustered on individual) reported in brackets. +=.1, *=.05, **=.01, ***=0.001")
)
## Calculating marginal effects and diff-in-diff
Coefplot_OLS_robust_moderators <- model_marginal_effects(OLS_robust_moderators) %>% mutate(model = moderators_names)
Coefplot_OLS_robust_moderators_polbuff <- model_marginal_effects(OLS_robust_moderators_polbuff) %>% mutate(model = moderators_names)
## Outputting Coefficient Plots
Coefplot_moderators <- bind_rows(Coefplot_OLS_robust_moderators_polbuff,
Coefplot_OLS_robust_moderators) %>%
pivot_longer(c(state_secretaries_NSDID, bureaucratic_buffer), values_drop_na = T) %>%
mutate(model = factor(model, levels = rev(str_replace_all(c("Baseline",
"government_PMparty_share_posts_end_of_year",
"years_since_DEP_start", "dep_end",
"multiple_L1", "age_group_year2"), output_names))))
print_plot(marginal_effect_plot(Coefplot_moderators, models = T), name = "Coefplot_moderators")
savehistory("C:/Users/bjormf/UiO Dropbox/Bjørn Forum/Apps/Overleaf/bufferingpaper_BJPolS/.Rhistory")
